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ABSTRACT 

The initial value problem for the first order linear wave equation in one space dimension is treated for two cases 
with specified initial data and grid, and data from solutions at t = 400 and t = 800 are presented, as prescribed for 
Problem 1 in Category 1 . Results are shown from computations with a sequence of recently developed high order and 
high resolution methods which combine Hermite interpolation, Cauchy-Kowaleskya recursion for time derivatives, 
and Taylor series time advancement [3, 4, 5], These methods have the same order of accuracy in time as in space. 
Results are shown from methods that range from third to nineteenth order. The stated problems with the prescribed 
coarse grid can be simulated with errors that are at the level of machine accuracy if the method is sufficiently high 
order. In addition, the growth of the maximum absolute error out to t = 100, 000 is given for simulations with the 
stated problem data. 


PROBLEM STATEMENT 


The initial value problem is for the first order linear wave equation in one space dimension. 


du 

dt 


du 

dx 


= 0, 


( 1 ) 


with initial data 

u(x,0) = (2 + cos[aa;])exp[-ln[2](^) 2 ]. (2) 

The two cases that are required are with a = 1.7 and a = 4.6, with solution data at t = 400 and t = 800. A uniform 
grid is supposed is to be used, with mesh size Ax = 1. Any method that captures the physical dynamics for this 
problem should be able to produce machine accuracy for time step size At = Ax, since this choice simply propagates 
the solution along the characteristic from one grid point to the next over one time step. The computations that are 
shown here use At = | Ax, so that the grid data is not simply propagated from one point to the next. 


SOLUTION METHOD 


Results are presented from a sequence of recently developed high order and high resolution methods which com- 
bine Hermite interpolation, Cauchy-Kowaleskya recursion for time derivatives, and Taylor series time advancement 
[3, 4, 5], These methods can be thought of in a variaety of ways. Consider an expansion of the solution in both space 
and time about the grid point (x*, tj), with 


Uij(x,t) 


l{x Xi^t tj ) — ^ ^ ^ 


a, b—0 
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in local coordinates, where the grid index (i, j) has been suppressed in the series expansion on the right. If this 
expansion is a solution of the governing equation, then the expansion coefficients must satisfy 


^a,b + 1 4“ Ua+l,b 0, Or tta,fc+l t/ a _ 

which reflects the recursion 

Q a +b+ l u ga+b+l u 

dx a dt b+1 dx a+1 dt b 

from the governing equation. Consequently, the solution for the model problem can be written in the form 
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which is just the method of characteristics. Note that 


^ ^ j Uq,x 
O !=0 

is an expansion in space of the known solution at the current time, and that 

OO 1 

Ui,j(0,t) = ^2 —ii a (-t) a 

a—0 


(7) 


( 8 ) 


is an expansion of the solution in time at the stencil center. These two local expansions in space and in time are 
connected by the recursion relations that are derived from the governing equations. The spatial expansion is just initial 
data on a noncharacteristic surface that is to be propagated by the governing equations, and the recursion relations 
are just the Cauchy-Kowaleskya recursions [2] from the governing equations. The key to applying this method to 
more complex systems is to obtain the recursion relationship for the local expansion coefficients from the governing 
equations, just as above in this particularly simple case. 

Particular realisations of this method truncate the local expansion to a given order for spatial interpolation, and 
then the expansion in time is used to advance the solution to the same order at the stencil center. We use Hermite 
spatial interpolation with two point stencils on a staggered grid system. Two uniform Cartesian grids are used, where 
each is offset from the other by half a grid step. At each grid point at the current time level tj we keep data for the 
dependant variable and all of the spatial derivatives up to order to, or numerically computed approximations for 


{ 


d k u 

dx k 


0 , . . . , to }. 
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On a two point stencil this provides data for an interpolant that is of order 2 to + 1, with consistant estimates at the 
stencil center for the derivatives of u up to order 2m +1. In local coordinates about the center of the two point stencil 

\{xj + x j+1 ) = x j+ i we have 


2m+l ^ 

Ui j(x, 0) « —u a x a , where u a 
' a\ 

oc—0 


d a u 
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The time derivatives for the solution and for its spatial derivatives up to order to are obtained from the space derivatives 
by the Cauchy-Kowaleskya recursion relationships for the governing equation. For the workshop problem. 


d a+(3 u , . , a d a+/3 u . , 


(-1 ) 0 U a+f} . 
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The local expansion in time can now be used to produce a Taylor series time expansion for the solution and its first to 
derivatives. In the case of the simple workshop problem, the time expansions are just 


d k Ujj 

dx k 


(0,A t) 


2m-\-l—k .j 

V —u a+k {-At) a , 

z — ' a! 

a— 0 


( 12 ) 


for k = 0, 1, . . . , to. Note that a derivative of order k is advanced in time only to order 2 to + 1 — k, but that the 
solution itself (. k = 0) is advanced to order 2 to + 1. Since the expansion center is halfway between two grid points 
on one of the two staggered grids, the time advanced solution data is provided for a grid point on the other offset grid. 
For these computations Taylor series time advance is with At = 2 Ax. Further details are in [3, 4, 5], 

A summary of the algorithm realisations that are used for the workshop calculations is presented in Table 1. The 
first column on the left presents the method name. The second column gives the order to of the derivatives that are 
kept at each grid point, with derivative order m = 1, 3, 5, 7, 9. The third column gives the order 2 to+ 1 of the method, 
from 3 rd to 19 , in both space and time. The fourth column gives the FLOP count for the number of multiplications 
that are required for each grid point for each time step. The higher order methods are clearly more complex, the 
number of FLOPS per grid point per time step for the 19 th order c2o9 method is 26.5 times the count for the 3 rd 
order c2ol method. On a fixed grid the higher order methods will always be more expensive to use. However, if it 
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method 

derivative order 

method order 

FLOPS/g.p/t.s. 

c2ol 

1 

3 

26 

c2o3 

3 

7 

108 

c2o5 

5 

11 

246 

c2o7 

7 

15 

440 

c2o9 

9 

19 

690 


Table 1: Summary of Methods Used for Workshop Results 


is necessary to stay below a fixed level of error, and if grid refinement and coarsening can be done, then the higher 
order methods will be able to produce a required accuracy with significantly coarser grids and fewer time steps than 
the lower order methods. The result is that the total FLOP count to reach a required level of accuracy for an entire 
simulation is significantly lower with high order methods than with lower order methods. In general, computational 
efficiency increases with order (see [5] for more details). 

Because of their structure, these high-order high-resolution methods can be called Hermite/Cauchy-Kowaleskya/Taylor 
(HCKT) methods. Because of their development from local series expansions interrelated by the recursion relation- 
ships for the governing equations, they have also been called Modified Expansion Solution Approximation (MESA) 
methods [1], HCKT methods are stable for CFL < 1, [8], These methods can be of very high order accuracy, 
with very high resolution, and the high order methods are typically very computationally efficient, requiring orders 
of magnitude less total FLOPS to obtain a predetermined error level when compared to standard methods such as 
compact differencing with Runge-Kutta (see [5]). We have successfully used this approach with the linearized Euler, 
heat, Klein-Gordon, Burgers and Navier-Stokes Equations. In all cases the systems can be written with only first order 
time derivatives. For linear constant coefficient systems in Cartesian coordinates, the HCKT methods can have exact 
propagators in time even in multiple space dimensions. In the case of hyperbolic systems, if the propagation in time 
is exact, so that the dynamic evolution of the local initial data is correct, then information is correctly propagated 
along characteristic surfaces. In this sense, these methods can be viewed as a correct generalization to multiple space 
dimensions of the method of characteristics in one space dimension. The HCKT methods have been used succesfully 
with high order accuracy at both radiation [9, 7, 10] and surface [6] boundaries. The HCKT methods have also been 
used on randomized grids with no loss of accuracy and efficiency. 

SUMMARY OF ERROR DATA FOR THE WORKSHOP COMPUTATIONS 


method 

order 

a = 1.7 
t = 400 

t = 800 

a = 4.6 
t = 400 

t = 800 

c2ol 

3 

1.00 

1.01 

1.00 

1.01 

c2o3 

7 

7.12 x 10" 3 

1.42 x 10 -2 

9.96 x lO” 1 

9.96 x 10" 1 

c2o5 

11 

2.35 x 10" 7 

4.69 x 10~ 7 

2.82 x 10 -2 

5.56 x 10" 2 

c2o7 

15 

3.26 x 10" 12 

6.73 x 10” 12 

8.39 x 10 -6 

1.68 x 10" 5 

c2o9 

19 

1.20 x 10 -12 

1.73 x 10" 12 

1.49 x 10" 9 

2.97 x 10" 9 


Table 2: Maximum Absolute Error for Workshop Problems From All Methods 


The maximum absolute error for the workshop computations is presented in Table 2. Recall that all of these 
computations are with Ax = 1 and At = The numerical domain is for —1000 < x < 1000. The first column 
presents the method name for the data in each row. The second column gives the order of the method in both space 
and time. The third through sixth columns give maximum absolute errors in the solution u: columns three and four 
with a = 1.7, and columns five and six with a = 4.6; while columns three and five give the error at t = 400, and 
columns four and six at t = 800. To the scale of a solution plot, the errors are invisible if they are O[10~ 3 ]. For 
this error tolerance, the seventh order c2o3 method provides an adequate solution for a = 1.7, but the eleventh order 
c2o5 method is still not quite powerful enough for a - 4.6. Note that for a = 1.7, the fifteenth and nineteenth 
order methods produce errors that are in the machine roundoff range. For a = 4.6, the ninteenth order c2o9 method 
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produces errors that are 0[ 10 9 ], Clearly, powerful methods that are high order and high resolution can produce 
accurate solutions for this problem. 

ADDITIONAL COMPUTATIONS ON A PERIODIC DOMAIN 

The maximum absolute error at t 800 as a function of grid size with a = 1.7 is presented in Table 3, and 
with a = 4.6 in Table 4. This data is from computations on the domain —100 < x < 100 with periodic boundary 
conditions. A periodic domain was chosen to avoid the need for a very large grid that would arise if the domain 
— 1000 < x < 1000 were to be used with a small mesh size. The data in these Tables does not confirm the order of 
accuracy of the fiftenth and nineteenth order methods, because both methods produce very small errors on the coarest 
grid that has been used (Ax = 1). The other methods have their orders of accuracy confirmed by this data. 


1/Ax 

c2ol 

c2o3 

c2o5 

c2o7 

c2o9 

1 

1.01 x 10+° 

1.42 x 10 -2 

4.69 x 10” 7 

6.73 x 10 -12 

1.73 x 10" 12 

2 

9.98 x lO’ 1 

1.25 x 10” 4 

2.19 x 10" 10 

1.01 x 10" 12 


4 

5.80 x 10" 1 

1.00 x 10" 6 

2.97 x 10" 13 



8 

1.07 x 10” 1 

7.89 x 10 -9 

3.46 x 10" 13 



16 

1.41 x 10 -2 

6.17 x 10" 11 




32 

1.77 x 10" 3 

7.21 x 10" 13 




64 

2.22 x 10" 4 






Table 3: Maximum Absolute Error by Grid Density at t = 800 for a = 1.7 With Each Method 


1/Ax 

c2ol 

c2o3 

c2o5 

c2o7 

c2o9 

1 

1.01 x 10+° 

9.96 x 10“ 4 

5.56 x 10" 2 

1.68 x 10 -5 

7.48 x 10" 10 

2 

1.00 x 10+° 

2.29 x 10" 1 

2.85 x 10" 5 

8.95 x 10" 10 


4 

1.00 x 10+° 

2.43 x 10 -3 

1.29 x 10 -8 

2.52 x 10~ 12 


8 

9.97 x 10" 1 

2.00 x 10” 5 

6.03 x 10- 12 



16 

5.21 x 10 -2 

1.58 x 10" 7 




32 

8.85 x 10 -2 

1.24 x 10" 9 




64 

1.15 x 10" 2 

9.72 x 10 -12 





Table 4: Maximum Absolute Error by Grid Density at t = 800 for a = 4.6 With Each Method 


1/Ax 

c2ol 

c2o3 

c2o5 

c2o7 

c2o9 

1 

1.66 x 10 7 

6.91 x 10 7 

1.57 x 10 8 

2.82 x 10 8 

4.42 x 10 8 

2 

6.66 x 10 7 

2.76 x 10 8 

6.30 x 10 8 

1.13 x 10 9 


4 

2.66 x 10 8 

1.11 x 10 9 

2.52 x 10 9 

4.51 x 10 9 


8 

1.06 x 10 9 

4.42 x 10 9 

1.01 x 10 10 



16 

4.26 x 10 9 

1.77 x 10 10 




32 

1.70 x 10 10 

7.08 x 10 10 




64 

6.82 x 10 10 

2.83 x 10 11 





Table 5: Total FLOPS required to compute to t = 800 by Grid Density With Each Method 


Tabel 5 presents the total FLOPS required to do the computations that are summarised in Tables 3 and 4, by grid 
size for each method. Recall that there are two half time steps per full time step, and the coarsest grid has 200 grid 
points, so that the total FLOPS to compute a solution on the coarsest grid out to t = 800 is 200 x 2 x 1600 x F, where 
F is the FLOP counts per grid point per time step from Table 1 for each method. The grid ratio is kept at At = | Ax 
for all computations. The Total FLOP counts can be used with the data in Table 4 to make some useful comparisons. 
Note that an error of 0[ 10 -2 ] is obtained by the c2ol and c2o5 methods with grid densities of 1/Ax = 16 and 1, 
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respectively. Table 5 shows that the third order method requires approximately 27 times more computational effort 
than the ninth order method. Similarly, an error of O[10 -5 ] is obtained by the c2o3 and c2o7 methods with grid 
densities of L/Ax = 8 and 1, respectively. Table 5 shows that the seventh order method requires approximately 16 
times more computational effort than the eleventh order method. Finally, an error of O[10 -1 °] is obtained by the c2o3 
and c2o7 methods with grid densities of 1 / Ax = 32 and 2, respectively. Table 5 shows that the seventh order method 
requires approximately 63 times more computational effort than the eleventh order method. This data clearly shows 
a trend of greater efficiency for higher order methods, irrespective of the error level. We may also look at the level 
of accuracy produced by the different methods with a similar computational effort. Table 5 shows that four of the 
five methods have a total FLOP count entry that is O[10 9 ]. The c2ol, c2o3, c2o5, and c2o7 methods exert this effort 
with the grids 1/Ax = 8, 4, 4, and 2, respectively, with corresponding errors that are 0[1], O[10~ 3 ], O[10~ 8 ] and 
O[10 _1 °], respectively. This crude comparison shows a truly startling range of accuracy for the same computational 
effort, from errors that are the size of the solution itself to close to the machine roundoff level, for methods from third 
to fifteenth order accurate in both space and time. If a fuller range of grid resolutions were used, both finer grids for 
the low order methods and coarser grids for the high order methods, then a more systematic study of computational 
efficiency could be done, as in [5]. Nevertheless, it is clear from even this limited data that the higher order methods 
reach a given error limit more efficiently than lower order methods, and that a given level of effort produces much 
more accurate results with higher order methods than with lower order methods. 


t 

a = 0 

a = 2 

a = 4 

a = 6 

a = 8 

a = 10 

1 x 10 3 

1.69 x 10" 12 

2.61 x 10” 12 

2.44 x 10 _1 ° 

6.15 x 10- 7 

1.16 x 10 -4 

6.39 x 10 -3 

1 x 10 4 

9.64 x 10” 12 

9.95 x 10~ 12 

2.44 x 10~ 9 

6.13 x 10 -6 

1.14 x 10" 3 

5.36 x 10 -2 

2 x 10 4 

1.06 x 10 -11 

1.30 x 10” 11 

4.86 x 10 -9 

1.23 x 10" 5 

2.28 x 10" 3 

1.03 x 10” 1 

3 x 10 4 

1.40 x 10" 11 

1.42 x 10 -11 

7.29 x 10" 9 

1.84 x 10" 5 

3.41 x 10" 3 

1.51 x 10” 1 

4 x 10 4 

1.97 x 10“ n 

1.34 x 10" 11 

9.72 x 10“ 9 

2.45 x 10" 5 

4.55 x 10 -3 

1.96 x 10" 1 

5 x 10 4 

2.66 x 10 -11 

1.76 x 10" 11 

1.21 x 10" 8 

3.06 x 10" 5 

5.68 x 10 -3 

2.38 x 10" 1 

6 x 10 4 

3.31 x 10" 11 

1.93 x 10" 11 

1.46 x 10 -8 

3.68 x 10 -5 

6.82 x 10 -3 

2.78 x 10" 1 

7 x 10 4 

2.42 x 10 -11 

1.70 x 10" 11 

1.70 x 10" 8 

4.29 x 10" 5 

7.95 x 10 -3 

3.17 x 10" 1 

8 x 10 4 

2.83 x 10 -11 

2.88 x 10 -11 

1.94 x 10 -8 

4.90 x 10" 5 

9.08 x 10“ 3 

3.53 x 10" 1 

9 x 10 4 

3.22 x 10" 11 

2.54 x lO^ 11 

2.19 x 10" 8 

5.51 x 10 -5 

1.02 x 10 -2 

3.87 x 10- 1 

1 x 10 5 

3.37 x 10 -11 

2.92 x 10- 11 

2.43 x 10" 8 

6.13 x 10" 5 

1.13 x 10 -2 

4.20 x 10" 1 


Table 6: Maximum Absolute Error Over Time for Variaous a With Ax = 1 for the c2o9 Method. 


Table 6 presents maximum absolute error date for computations with the nineteenth order c2o9 method. These 
computations are all on the domain —100 < x < 100 with periodic boundaries, and on the coarsest grid with Ax = 1 
and At = \ . The compuatations are for a range of a from 0 to 10. Note that the computations run out to t = 100, 000, 
or 125 time further than for the workshop problem. The data shown is for the absolute maximum cumulative errors 
for times from t = 1000 to t = 100, 000. Note that the errors for a = 0 and a = 2 are very close, starting in the 
machine roundoff range at t = 1000, and increasing very slowly out to the final time t = 100, 000. The errors for 
the other four values of a increase much more rapidly, by the same order as the increase in time from the first to last 
entry. These computations give a more complete sample of the accuracy of the nineteenth order c2o9 method. Note 
that for a < 4, the error growth can be projected to later times, and suggest the possibility of a scarcely visible error 
that is O[10~ 3 ] at t = O[10 10 ]. This data shows that these high order and high resolution methods can be used to 
accurately propagate relatively small scale data over large distances or long times, or for computations where there is 
a large difference between the scale of the problem as a whole and the scale of the smallest information that must be 
resolved and accurately propagated. 


CONCLUSIONS 

New high order and high resolution methods have been applied to a benchmark problem for the fourth CAA 
workshop. These HCKT of MESA methods combine Hermite differencing with high order derivative data on two 
point stencils, Cauchy-Kowaleskaya recursion for computing time derivatives from space derivatives, and high order 
Taylor series time advancement between staggered grids. These methods all have the same order of accuracy in space 
and in time. Methods from third to nineteenth order were used to produce the results in this paper. The sample 
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computations show that errors down to the machine roundoff level can be obtained by a sufficiently accurate HCKT 
method even on a coarse grid such as has been specified for the workshop problem. Data from additional computations 
show that the more accurate methods can reach a stated error limit with much greater computational efficiency than 
lower order methods, or that they can maintain a much lower error for a given level of computational effort than 
lower order methods. Sample calculations show that the higher order methods can produce very accurate results for 
extremely long times. 
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